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static form for one dimensional PDEs 


Maryam Arab Ameri 


Abstract 


In this paper, we propose an adaptive mesh approach for time dependent parial 
differential equations, based on a so-called moving mesh PDE(MMPDE) and level 
set method. It means that the velocity of mesh nodes is calculated by MMPDE 
and is employed as veocity in the level set equation. Then, at each time level, the 
mesh points are considered as the level contours of the level set function. Finally 
the method is merged with local time step technique. 
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1 Introduction 


In this paper, we discuss a class of adaptive mesh algorithms for solving time- 
dependent partial differential equations(PDEs) whose solutions have large 
variations over a given physical domain, such as shock waves, boundary layer 
and interior layer. 

This method is based on the level set concept. the level set methods are 
powerful numerical techniques for tracking the evolution of interfaces moving 
in complex ways. The level set methods were used for the first time to repre- 
sent the evolution of surfaces implicity by Osher and Sethian [11]. Due to its 
many advantages, this approach has been used for many cases [1, 5, 9, 13, 19], 
for example, it has been used for mesh generation around a body(inside or 
outside) by the level set function [12]. Liao et al. presented some points about 
using the level set functions for moving grid based on deformation method 
[6]. In this paper, we also formulate an adaptive mesh method which is based 
on the level set method. This method is combined with a class of moving 
mesh methods which employs a moving mesh partial differential equation, 
MMPDE, to perform mesh adaptation. The MMPDE is formulated in terms 
of coordinate transformation or mapping,[3, 4, 8, 16, 17]. Several moving 
mesh equations(MMPDEs) based upon equidistribution principle have been 
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derived in [3], where in Section 2, we briefly review some of them. Finally 
after introducing the new method, we present local time stepping refinement 
technique based on the slope of the level set function for increasing the effi- 
ciency of the moving mesh method. The mentioned method is used to solve 
scalar field satisfying 

u(x,t) = L(u), (1) 


where L is a differential operator defined on physical domain 2, a < x < b 
and 0<t< Ty. 


2 Moving Mesh Methods 


In this section, we develop a new moving mesh method based on the level set 
concept. We start with a review of the equidistribution principle(EP)|[3, 15] 
and derive MMPDEs from that principle in subsection 2.1,then a description 
of adaptive grid based on level set method is given in the subsection 2.2. 


2.1 MMPDE 


The evolution of moving computational grid can be viewed as a discretiza- 
tion of a one-to-one time dependent coordinate mapping. Let « and € denote 
the physical and computational coordinates, with domains 2 and §2., respec- 
tively. Without loss of generality, both of them are assumed to be in [0, 1]. 
Define a coordinate transformation by: 


a = x(€,t), €€ [0,1], 2(0,t) =0, x(1,t)=1. (2) 


The computational coordinates is discretized on a uniform mesh given by 


&= 7=O0, 1,000 ,N; 


N ) 
where N is a certain positive integer and the corresponding mesh in x denoted 
by 

0 = 20 < a(t) < to(t) <-++ << ay_i(t) << ay = 1. 


A major factor of moving mesh approach is the monitor function, p(2,t), 
which is chosen to be somehow a measure of the solution error. For a given 
monitor function, the mesh point locations, x; (t), could be required to satisfy 
the following equidistribution principle (EP) for all the values of time t¢ [3]: 


x; (t) 1 1 ii 
| px, thd = — i p(x, td = —0(t), 


j—1(t) 
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or equivalently 
x5(t) j 
[ ole.tiae = Lote) = 6010) (3) 
0 


By differentiating (3) with respect to € once and twice,, we obtain two dif- 
ferential forms of EP, 


() 
pla(€,t)) seals, t) = O(t), (4) 
3) 3) 
pg (PE), Heats, 6) = 0. (5) 


These EPs,(3),(4) and (5), which do not contain the node speed i(€,t),are 
called quasi-static EPs(QSEPs). Related to these QSEPs, various MMPDEs 
were derived in [3] by taking the mesh to satisfy the above EP (Eq. (5)) at a 
later time t+ 7 instead of at t. In this case the mesh should satisfy 

e t t a t =0 6 
perires, +7), +r) ages, +7) _ ( ) 
where the parameter 7 is called a relaxation time. By expanding the term 
Za(€ ,t+7) in Taylor series and dropping certain higher order terms, various 
MMPDEs can be obtained. Two of them which are used in our work are 
MMPDE5 and MMPDE6: 


. 10, 02 

MMPDES: —i=—=7.(p52); (7) 
Ora 10, Ox 

MMPDE6: 55 = ~~ 5¢(P5p): (8) 


These MMPDEs and also the rest of them not only force the mesh ,x(€,t), 
toward equidistribution principle but also prevent the mesh from crossing. 
More specifically, the term —(4)( 8 (03) plays the fundamental role of a 
correction term to make the mesh equidistribute the monitor function and it 
is stabilizing term for the mesh trajectories. 

The monitor function, p, is chosen such that the mesh points are concentrated 
in regions where more accuracy is needed, and so p is usually taken to be 
some measure of the solution error estimated from discrete solution values. A 
commonly used monitor function is p = \/1+ au? (a is regularizing factor), 
which equidistributes the arclength of the solution u. 

In this work, we discretize MMPDE5 and MMPDEG6 by centered finite dif- 
ference in space, which yields for MMPDE5 


t=, j=0,1,...,N 


and for MMPDE6 
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Bi+1— 285 +2j1=—-—, 9 =O, Djoe cy NV: 


The quantity E; represents a centered approximation to the term on the right 
hand side of MMPDE5 and MMPDE6 given by 


E; _ Myj+1/2(@j41 — xj) _ Mj-1/2(2; = oi), a = 0, 1, at .N 


where Mj41/2 = 5(Mj + Mj41),Mj; = M(u;), and u; © u(a;,t) is an ap- 
proximation of the solution at grid point z;. 


2.2 Moving Mesh Based on Level Set Approach 


In this section, at first a new moving grid method is formulated, and then 
the algorithm of this method is presented. 

Essential point in all applications of the level set method is using the implicit 
representation. This point is also used in grid generation in such a way that 
the mesh nodes are represented implicitly by the level contours of the level 
set function. In this method, at each time level we construct the level set 
function,y(z,t), and the mesh points z;,j =0,1,...,.N, are obtained by the 
level contours of w(a,t), that means the mesh points are: 


{xz; | o(a,t) =c;, 7 =1,2,...,N} (9) 


where c; is j-th component of a constant vector C = (0, a <, ...,1) and N 
is the number of mesh points. 

In fact, in each time, we seek a level set function, € = w(a,t): 23 2Q 
, which maps x; to c; = d, j = 0,1,...,N and satisfies the well-known 
equidistribution principle . 


Jy = —— 
P| Qe’ 


(10) 
where J is the Jacobian of the coordinate transformation x = w!(€,t) 
2. + 2, and p = p(x,t) is a given monitor function, and o = a(t) 
Jo ela, thdz. 

For updating the level set function or equivalently for updating the position 
of mesh nodes, the well-known level set equation is applied, 


U(x, t) + v_(x,t) = 0. (11) 


The ideal initial condition for the above equation is ~(#,t) = x, because in 
the most mesh adaptation algorithms, the purpose is to convert a uniform 
mesh at initial time to an equidistributed mesh at the next time levels. The 
mentioned initial condition gives a uniform mesh at t = 0. Also the following 
boundary conditions are considered for the level set equation, 
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W(0,t) = 0, ¥(1, t) = 1. 
As we mentioned at each time level, any 
Wet) =a (12) 


has a correspondence point on x-axis which is one of the new mesh nodes. 
By differentiating (12) with respect to t, we get 


Yr (x,t) + bo(a, t)t = 0, (13) 
and when compared with (11) and (13), we get 
£=0, (14) 


that means, the nodes velocity is equal to v. 

So for calculating the nodes velocity in the level set equation, the MMPDE5 
or MMPDE6 should be applied which generates an equidistributed mesh 
through this algorithm. 

At each time level, after determining new nodes, the solution of PDE should 
be determined. For this purpose, let 


U(w(a,t),t) = u(a,t), (15) 
then 7 
U; => Une + Ut 


where wu; = L(u) by (1) and & = v by (14). The derivatives in L(w), such as 
Ux, Ux are also transformed. For example, from (15), we have 


- U, U, 2 ~ 
Uy Ur®y Ue y = pl a a, 
Th Qe] o 


Coy —Ualvy _ (eee 

(ty)? o 
The higher derivatives can be derived similarly. The transformed equation 
for U(w,t) takes the form of 


)? (Opp — Ua typ). 


Taps = en ConA +Urlyy > Ure = 


where L is a differential operator in 7. Finally (16) will be solved on a uniform 
grid. 
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2.3 Algorithm of Adaptive Level Set Method(ALSM) 


In this section, we provide an algorithm to solve the given PDE (1) in a 
moving grid based on the level set function. 

Step1: Enter the initial time, to, the final time, 7, the length of time step, 
At, and the number of mesh points, N. 

Step2: Set i= 0. 

Step3: If tp = 0, set 


a 
KO = ig = ae =0,1,...,.N}, 
pio) = xo), 
ytto) = yfto) = u(X) ¢9), 
ito) = py(to) 


Step4: Determine p(u(z,t;)) by the solution u being calculated. 
Step5: Compute mesh velocity, v = «, either by MMPDE5: 


i=, j=0,1,...,N 
or by MMPDE6: 
E; 
jy — 2bj + ja = ——. j=0,1,...,N 


Step6: Update the level set function by the following equation: 
U(a, t) =e Wx (x, t).v = 0, 
which by forward finite difference discretization in time converts to: 


W(tit1) = Uti) + At(vdes)e:, 


with the boundary conditions w(t; + 1) =0 and qn (ti +1) = 1. 

Step7: Define the inverse of w(a,t) for updating the new nodes, X(t;+1), in 
current time. 

Step8: Determine the values of u on the current time level by the described 
procedure in previous subsection and solve (16) or 


O(tiz1) = U(t;) + At.L(U (ti), 


Step9: Set i=i+1. 
Step10: If t; < 7 go to Step3, else break. 
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3 Local Time Step Refinement 


Local time stepping for one-dimensional conservation laws was first proposed 
by Osher and Sanders [10]. Tan et. al. proposed variable time stepping for one 
and two dimension for conservation laws, which uses semi-dynamic moving 
mesh method such that the CFL condition is used to obtain time interval 
refinement to compute the solution [18]. Also, Soheili and Salahshour used 
this approach for the blow-up problems [14]. 

In this section, we also present the details of the local time stepping refine- 
ment(LTSR) which is performed based on the level set function. 

Let initial time steps of the problem have the form 


to =0 > t) = tpt At > ... > th = to NAF > thai = to+(n4+1)At > ... > Ty, 


where At is a specified value for the time step and is constant through using 
this procedure. On the first interval [to,t1], set Ato = a ae where kg € N is 
constant and depends on the slope of the level set function, (the method of 
determining ko will be described). We have 


to4(ko—i) Ato =to+ (ko = i) Ato, i ko, ko = i. wae ,0 


so the time integration on [to,t1] involves ko sub-steps such that 


to =0 >to + Ato >... 3 to + ko Ato = fh. 


Generally, suppose that we are at time level t = t, and we want to move 
towards t = ty41. Similarly consider At, = fnta—tn such that 


tnt (kn—t)Atn =tyht+ (Kin _ i)Atn, t= kn, kin — 1, fey 0 


that means, on interval [tp,tnr41], we have 


th 3 th + Aty >... > th + ky Aty = tris. 


This process continues up to Tr. Now, ko, ki,...,kn,... are determined by 
the following process. At the first step, we start integrate the problem from 
to to ty = to + At without any LTSR, we call this process prediction step, 
then we determine the level set function, the new mesh nodes of adaptive 
mesh and solution of PDE at t,. According to the property of the mentioned 
adaptive mesh method in previous section, ”for larger slope of the level set 
function, more mesh nodes are concentrated”. Thus according to slope of 
the level set function, we define the interval [t,t,]. So the slope of pi) 
on [%;,%;41], 7 = 0,1,...,N is calculated, then for having more efficient 
solution, [to,t1] is subdivided to ko parts. Judicious choice for kg can be 


ko = min{max[slope(\”), 10)}}, 7 =0,1,...,N, 
j 
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and the time integration is done on this sub-interval again but with Atp = 
? After obtaining the solution at t, we act on [t,,t2] similarly and de- 
termine the solution at to. 

Generally, after calculating the solution at t,, at first we have a prediction 
process concluding time integration from t, to tn4+1 with pre-determined 
PP) os (apy) d= 
0,1,...,N and subsequently k, = min{max[slope(1s"*”)], 10}. Then the 


value of At as time step for knowing the slope of wh 


time integration is performed on this sub-interval k,, times with local time 
step At,,. 


3.1 Algorithm of Adaptive Level Set Method with LTSR 


Step1: Enter the initial time, to, the final time, 7’, the length of time step, 
At, and the number of mesh points, N. 

Step2: Set i= 0. 

Step3: If tp = 0, set 


pfto) = pylto) = a (to) u( Xo), to). 


Step4: Run ALSM one time with to = ¢;,T, = t;4, and At, then calculate 

pitti) X tits) and Ui+1), 

Step5: Calculate the slope of p+) = p+) on [zj,a;4:], jg =0,1,...,N 
(i+1) 


and subsequently k; = min{max(slope(y; " ’)], 10}. 
j 


Step6: Set At; = Ata 

Step7: Run ALSM kj times with to = tj, Ty = tiz1, and At = At;, then 
again calculate i+), X(+1) and Ui), 

Step8: Set i=i+1. 

Step9: If t; < Ty go to step4, else break. 


4 Numerical Experiment 


In this section, we implement the present work for two numerical examples 
where u(x,t) is a given function. We use the arclength monitor function. In 
order to obtain an accurate and non-oscillatory solution, it is necessary to 
smooth the mesh points. Following [5], we have applied smoothed monitor 
function as below: 
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iti ki 
Det a ont) i 


Pi ~~ ita 3 
ae er eat 


? 


where ip is a nonnegative integer and 7¥ is a positive constant. In this paper 
we select the mentioned formula for smoothing monitor function with ip = 1 
and y = 1. We apply MMPDE5 and MMPDE.6 in the evolution equation of 
the level set function for obtaining nodes velocity. 


(a) (b) 


Fig. 1: The mesh trajectory and solution of the first example for 0 < t < 1.0 with adaptive 
level set method(MMPDE5) of 21 mesh point,(w = 0.5, A = 0.4,¢ = 0.01, 6 = 0). 


Example 4.1. Consider the Burger’s equation as first example, 
Ut + Uz — EUge = O,7 
where its exact solution is: 


1) et A+ (p—A)es @-Ht-B) 
tla, = 1+ e2 (w—ut-B) Y 


We select pp = 0.5, A = 0.4,¢ = 0.01 and 8 = 0. This problem is characterized 
by moving discontinuities(specially when ¢ is very small), that means the dis- 
continuities move in time, and so the solution at a particular point in space 
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(a) (b) 


Fig. 2: The mesh trajectory and solution of the first example for 0 < t < 1.0 with adaptive 
level set method(MMPDE6) of 21 mesh point,(w = 0.5, A = 0.4,¢ = 0.01, 6 = 0). 


Table 1: The error of ALSM(with MMPDES) for the first example at t = 1.0, obtained by 
the arclength monitor function with a = O(uniform mesh) and a = 2(moving mesh)without 
LTSR and with LTSR, (4 = 0.5,A = 0.4,¢ = 0.01 and 6 = 0). 


MMPDE Time stepping Type of mesh N L°™-error CPU time 
21 0.1924 0.504 
5 without LTSR fixed mesh 31 0.0782 0.612 
a=0) 41 0.0564 0.8807 


21 0.0711 0.7548 
5 without LTSR Moving mesh 31 0.0602 0.8315 
a= 2) 41 0.0587 0.9812 


21 0.0510 0.917 
5 with LTSR Moving mesh 31 0.0347 1.019 
a= 2) 41 0.0298 1.079 


can change very rapidly. The solution of such a problem on a fixed uniform 
spatial mesh, needs very small time step in order to have sufficient accuracy, 
but using the adaptive mesh for finding the solution of this problem improves 
both the accuracy and the efficiency. 
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solution at +0 solution at tt 


(a) (b) 


Fig. 3: The solution at t = 0 with uniform mesh and at t = 1 with adaptive level set 
method(MMPDE5) of 21 mesh point for the first example, (u = 0.5,A = 0.4,¢ = 0.01, 8 
0). 


Table 2: The error of ALSM(with MMPDE6) for the first example at t = 1.0, ob- 
tained by the arclength monitor function with a = 0 (uniform mesh) and a = 2(moving 
mesh)without LTSR and with LTSR, (u = 0.5,A = 0.4,e¢ = 0.01 and 6 = 0). 


MMPDE Time stepping Type of mesh N L%-error CPU time 
21 0.0645 0.8815 
6 without LTSR fixed mesh 31 0.0578 0.9205 
(a = 0) 41 0.0512 0.9876 


21 0.0441 0.031 
6 without LTSR Moving mesh 31 0.0315 1.108 
(a = 2) 41 0.0274 1.230 


For this reason, we have used new adaptive mesh for this problem and 
in order to demonstrate the efficiency of the adaptive level set method 
(ALSM)(which has combined with MMPDE5S), and also to compare ALSM 
with or without LTSR, some results are presented in Tablel. These results 
show L°°-error and CPU times for different number of mesh points. These re- 
sults certify higher accuracy of the mentioned method(ALSM) in comparison 
with using uniform mesh. Besides, it shows that using the adaptive method 
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solution at +0 solution at +1 


(a) (b) 


Fig. 4: The solution at ¢ = 0 with uniform mesh and at t = 1 with adaptive level set 
method(MMPDE6) of 21 mesh point for the first example, (u = 0.5,A = 0.4,¢ = 0.01, 8 
0). 


with LTSR gives better results. Similar information, are given in Table 2 
when the MMPDE6 is applied for calculating nodes velocity. Figure 1(a) 
shows mesh trajectories for 0 < x < 2,0 < t < 1.0 which has been derived by 
using the new method with MMPDE5, also the solution of PDE, u, has been 
plotted on the moving mesh for 0 < t < 1.0 in Figure 1(b). The spatial do- 
main is divide into 21 mesh points. Like above, mesh trajectories and solution 
of Example 4.1 have been plotted by using the new method with MMPDE6 
in Figure 2. Also in Figure 3 we plotted the solution of PDE at t = 0 on 
a uniform initial mesh and at ¢ = 1 on a moving mesh with MMPDE5 and 
similarly, the solution at t = 0 and t = 1 has been plotted with MMPDE6 in 
Figure 4. 

For demonstrating the efficiency of our method we compare this method with 
another moving mesh method. For this purpose, we consider moving element 
free Petrov-Galerkin viscous method, MEFP-GVP, where introduced in [2], 
for comparing under equal conditions we consider the parameters ju, A,¢ and 
@ in the Burger’s equation according to [2], (u = 0.5,A = 0.4,¢ = 1/15 and 
(G = 0.16). Also, we solve this equation for x € [0,1] and t € [0,0.51] again. 
With above conditions and with N = 7 grid points the L°°-error is 0.03 by 
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(a) (b) 


Fig. 5: The mesh trajectory and solution of the first example for 0 < t < 1.0 with 
adaptive level set method(MMPDE5) of 21 mesh point,(j = 0.5, = 0.4,e = 1/15, 8 = 0). 


MEFP-GVM but in ALSM we have L°- error= 0.01. This result shows the 
preference of our method. In addition, the CPU time in our method is very 
smaller than the MEFP-GVM. Also we plotted grid motion and the solution 
for new conditions for 21 grid points in Figure 5. 


Table 3: The error of ALSM(with MMPDE5) for the second example at t = 1.0, obtained 
by the arclength monitor function (a = 2) without LTSR and with LTSR. 


MMPDE Time stepping N L-error CPU time 
21 0.0060 0.7213 
5 without LTSR 31 0.0052 0.8507 
41 0.0050 0.9759 


21 0.0051 0.906 
5 without LTSR 31 0.0032 1.015 
41 0.0029 1.076 


Example 4.2. Consider the equation 
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(b) 


Fig. 6: The mesh trajectory and solution of the second problem for 0 < x < 1 and 
0 <t< 1.0 with adaptive level set method(MMPDE5) of 21 mesh point. 


Table 4: The error of ALSM(with MMPDE6) for the second example at ¢ = 1.0, obtained 


by the arclength monitor function (a = 2) without LTSR and with LTSR. 


MMPDE Time stepping N 


L°-error CPU time 


21 0.0023 0.8706 
6 without LTSR 31 0.0015 0.9844 
41 0.0009 1.029 
21 0.0048 1.025 
6 without LTSR 31 0.0030 1.098 
41 0.0024 1.221 
Ut = Ura, 


with the exact solution: 


u(x,t) =e" ‘sin rx 


This problem has been used in [3] as a numerical example and also in [15] to 
study the moving mesh with variable relaxation time. 
Since u(x,t) — 0 in the limit as t + +0, then for typical arclength monitor 


function 
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(a) (b) 


Fig. 7: The mesh trajectory and solution of the second problem for 0 < x < 1 and 
0 <t< 1.0 with adaptive level set method(MMPDE6) of 21 mesh point. 


M=vV/1+4+au2 


we have M(x,t) > 1 as t ~ +00 and so the equidistributed mesh should 
tend to a uniform mesh in space. 

Similar to Tables 1 and 2, some results about the second example are given in 
Tables 3 and 4 which certify the quality of the new mesh adaptive algorithm. 
Figure 6(a) shows the mesh trajectories for the above problem using adaptive 
level set method with MMPDE5 and in Figure 6(b), the numerical solution 
of PDE is plotted for 0 < t < 1.0 by the new method with MMPDE6 have 
been plotted. We also plotted the solution at ¢ = 0 on uniform initial mesh 
and at the times t = 0.3,0.5 and t = 1 on moving mesh in Figure 8. As we 
expect, the adaptive mesh at t = 1 tends towards uniform mesh. 

In the second example, because of the smoothness of the solution there is not 
any major difference between the uniform and the adaptive mesh and also 
between the adaptive mesh with MMPDE5 and MMPDEG6. In this example, 
ALSM with LTSR does not yield any different result from ALSM without 
LTSR. 


Example 4.3. We use the following Burger’s equation as the third example, 


Ut + UU — EUger =O, OS & < xp. 
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solution at 0.3 


,0.3,0.5 and t = 1 with its corresponding 


Fig. 8: The solution of the second PDE at t = 0 


mesh nodes. 


wn 


Fig. 9: The mesh trajectory and the computed solution of the third example up to t = 1 
for 0 < x < 2 with adaptive level set method( MMPDE5) of 21 mesh points. 
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The initial condition is considered as follows [7]: 


2 
u(x, 0) = | sin(—*)), 
IR 


The boundary conditions are chosen in the form 
u(0,t) = u(xr,t) = 0. 


We select ¢ = 107° and rr = 2 in this example. This example has been 
solved by the presented method up to t = 1. 

The mesh trajectory has been plotted along with the nodes of moving grid 
with N = 21 mesh points in different times in Figure 9. Also the computed 
solution of the third example has been plotted in Figure 9. It is clear that the 
mesh moves correctly, because the solution has large variations about « = 0.5 
and x = 2 and the mesh also concentrates around these areas. 


5 Conclusion 


In this paper, we have developed a static moving mesh method based on the 
level set approach for solving one dimensional time dependent PDEs, such 
that for representing the nodes of adaptive mesh, the level set equation is 
used. Also, we proposed a strategy for local time step refinement where the 
local time step is selected by the level set function. This adaptive method is 
among the static moving mesh. We plan to investigate the dynamical moving 
mesh method based on the level set function. 
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